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Abstract. In this paper, we extend Stein’s method to products of 
independent beta, gamma, generalised gamma and mean zero normal 
random variables. In particular, we obtain Stein operators for mixed 
products of these distributions, which include the classical beta, gamma 
and normal Stein operators as special cases. These operators lead us to 
closed-form expressions involving the Meijer G-function for the proba¬ 
bility density function and characteristic function of the mixed product 
of independent beta, gamma and central normal random variables. 


1 Introduction 

In 1972, Stein [39] introduced a powerful method for deriving bounds for 
normal approximation. The method rests on the following characterisation 
of the normal distribution: W iV(0, a 2 ) if and only if 

E[A z f(W)} = 0 (1.1) 

for all real-valued absolutely continuous functions / such that E|/'(Z)| < oo 
for Z ~ N( 0,<x 2 ), where the operator Az , given by Azf{x) = u 2 f\x) — 
xf(x), is often referred to as the Stein operator (see, for example, Ley et al. 
[24]). This gives rise to the following inhomogeneous differential equation, 
known as the Stein equation: 

Azf(x) = h{x)-Eh(Z), (1.2) 

where Z ~ 1V(0, a 2 ), and the test function h is real-valued. For any bounded 
test function, a solution fy to (1.2) exists (see Stein [40]). Now, evaluating 
both sides at any random variable W and taking expectations gives 

E[A z fh(W)\ = Eh(W) - E h(Z). (1.3) 

Primary 60F05, 60E10 

Keywords and phrases. Stein’s method, normal distribution, beta distribution, gamma 
distribution, generalised gamma distribution, products of random variables distribution, 
Meijer G-function 

1 

imsart-bjps ver. 2011/11/15 file: ngb_bjps_2016.tex date: December 14, 2016 



2 


R. E. Gaunt 


Thus, the problem of bounding the quantity Kh(W) — E h{Z) reduces to 
solving (1.2) and bounding the left-hand side of (1.3). This is of interest be¬ 
cause there are a number of probability distances (for example, Kolmogorov 
and Wasserstein) of the form d%{£(W), C(Z)) = sup fte ^ \Eh(W) — Kh(Z)\. 
Hence 

d H (C(W),C(Z))< sup \E[Azf(W)}\, 
f G-F(ft) 

where T^K) = {fh ■ h G 77} is the collection of solutions to (1.2) for func¬ 
tions h £ This basic approach applies equally well to non-normal limits, 
although different Stein operators are needed for different limit distribu¬ 
tions. For a nice account of the general approach, we refer the reader to Ley 
et al. [24], 

Over the years, Stein’s method has been adapted to many other distribu¬ 
tions, such as the Poisson [4], exponential [3], [31], gamma [17] [25], [29] and 
beta [8], [19]. The first step in extending Stein’s method to a new probability 
distribution is to obtain a Stein equation. For the Beta(a, b ) distribution with 
density B ^ a fe ^ x a ~ 1 (l — x) b_1 , 0 < x < 1, where B(a, b) = r(a)r(6)/r(a + b ) 
is the beta function, a Stein operator commonly used in the literature (see 
[8], [19] and [36]) is 

-4beta fix) = x{l - x)f'(x) + (a - (a + b)x)f(x). (1.4) 

For the r(r, A) distribution with density j^yX r_1 e _Aa: , x > 0, the Stein 
operator 

Aamma f (x) = xf'(x) + (r - A x)f(x) (1.5) 

is often used in the literature (see [6] and [25]). In this paper, we extend 
Stein’s method to products of independent beta, gamma, generalised gamma 
and central normal random variables. In particular, we obtain natural gener¬ 
alisations of the operators (1.2), (1.4) and (1.5) to products of such random 
variables. 

1.1 Products of independent normal, beta and gamma random 
variables 

The theory of products of independent random variables is far less well- 
developed than that for sums of independent random variables, despite ap¬ 
pearing naturally in a various applications, such as the limits in a number of 
random graph and urn models (Hermann and Pfaffelhuber [22] and Pekoz et 
al. [32]). However, fundamental methods for the derivation of the probabil¬ 
ity density function of products of independent random variables have been 
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developed by Springer and Thompson [37]. Using the Mellin integral trans¬ 
form (as suggested by Epstein [11]), the authors obtained explicit formulas 
for products of independent Cauchy and mean-zero normal variables, and 
some special cases of beta variables. Building on this work, Springer and 
Thompson [38] showed that the p.d.f.s of the mixed product of mutually 
independent beta and gamma variables, and the products of independent 
central normal variables are Meijer G -functions (defined in Appendix B). 

The p.cl.f. of the product Z = Z\Z 2 ■ ■ ■ Zjy of independent normal random 
variables Z % r-*j N( 0, of), * = 1,2, ...,1V, is given by 


p(x) 


1 r N ,0 ( X 2 

(27 t) n / 2 <t °’ N \2 N a 2 



iel, 


( 1 . 6 ) 


where a = cr\02 ■ ■ ■ <Jn- If Z has density (1.6), we say Z has a product 
normal distribution, and write Z PN (N,(t 2 ). The density of the product 
X\ - ■ ■ X m Y 1 ■■■Y n , where X t r^i Beta (ai,bi) and Yj ~ T(rj, A) and the X t 
and Yj are mutually independent, is, for x > 0, given by 


p(x) = KG 


m+n, 0 
m,m+n 


where 



a\ + b\ — 1, a 2 + 62 — 1 ,..., a m + b m — 1 \ 
a l-l, a 2 - 1, ■ ■ • ,a m ~ 1, n - 1,..., r n - l) ' 

(1.7) 


K = X n l 


1=1 


T(oi + bi) 


n 


and we adopt the convention that the empty product is 1 . A random variable 
with density (1.7) is said to have a product beta-gamma distribution. If (1.7) 
holds with n = 0 , the random variable is said to have a product beta distri¬ 
bution, denoted by PB(ai, 61 ,..., a m , 6 m ); if (1.7) holds with m = 0, then 
we call this a product gamma distribution, denoted by PG(ri,..., r m , A). We 
also say that a product of mutually independent beta, gamma and central 
normal random variables has a product beta-gamma-normal distribution. 

For the product of two normals, (1.6) simplifies to 


p(x) = —-—x € M, 

7rcriC72 V (T l fT 2/ 


where Kq(x) is a modified Bessel function of the second kind (defined in 
Appendix B). For the product of two gammas, (1.7) also simplifies (see 
Malik [27]): 


p(x) 


2 \ri+r2 


r (ri+r 2 )/2-l tv- 

Us J-\.y —7*2 


(2\y/x), 


x > 0 . 
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Nadarajah and Kotz [28] also give a formula, in terms of the Kummer func¬ 
tion, for the density of the product of independent beta and gamma random 
variables. However, in general, for 3 or more (mixed) products of indepen¬ 
dent beta, gamma and central normal random variables there are no such 
simplifications. 

Pekoz et al. [33] extended Stein’s method to generalised gamma random 
variables, denoted by GG(r,X,q), having density 

p( x ) = Y^r- ) x r - 1 e-W 9 , x>0 . (1.8) 

For G ~ GG(r, X,q), we have EG k = X~ g T((r + k)/q)/T(r/q) and in par¬ 
ticular E G q = Special cases include GG(r, A, 1) = r(r, A) and also 
GG(1, (\/2 (t) _ 1 ,2) = HN(fi 2 ), where HN(cr 2 ) denotes a half-normal random 
variable: \Z\ where Z ~ IV(0, cr 2 ) (see Dobler [7] for Stein’s method for the 
half normal distribution). In this paper, we also extend Stein’s method to 
the product of generalised gamma random variables, GG(r*, A, q), denoted 
by PGG(n, ... ,r n ,X, q). 

1.2 Product distribution Stein operators 

Recently, Gaunt [13] extended Stein’s method to the product normal distri¬ 
bution, obtaining the following Stein operator for the PN(IV, cr 2 ) distribu¬ 
tion: 

■Azf(x) = cr 2 A N f(x ) - xf(x), (1.9) 

where the operator An is given by A^f{x) = x~ l T N f{x) and Tf(x ) = 
xf'(x). The Stein operator (1.9) is a IV-th order differential operator that 
generalises the normal Stein operator (1.2) in a natural manner to products. 
Such Stein operators are uncommon in the literature with the only other 
example being the IV-th order operators of Goldstein and Reinert [18], in¬ 
volving orthogonal polynomials, for the normal distribution. Very recently, 
Arras et al. [1] have obtained a IV-th order Stein operator for the distribu¬ 
tion of a linear combination of N independent gammas random variables, 
although their Fourier approach is very different to ours. Also, in recent 
years, second order operators involving /, f and f" have appeared in the 
literature for the Laplace [34], variance-gamma distributions [10], [14], gen¬ 
eralized hyperbolic distributions [15] and the PRR family of [32]. 

One of the main contributions of this paper is an extension of the product 
normal Stein operator (1.9) to mixed products of beta, gamma and normal 
random variables (see Propositions 2.3, 2.4 and 2.5). The Stein operators for 
these product distributions (given in Table 1) are higher order differential 
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operators, which cannot be readily be obtained via standard methods, such 
as the generator method of Barbour [2] and Gotze [20] and the density 
method of Stein et al. [41]. 

To obtain our product Stein equations, we use a conditioning argument to 
develop an algebra of Stein operators. In our proofs, we shall use the differ¬ 
ential operators T r f(x) = xf\x) + rf(x) and Bn,...,r„/(s) = T rn ... T ri f(x ) 
(note that Tq = T). It should be noted that whilst we restrict our attention 
to mixed products of betas, gammas and normals, we expect that the proofs 
techniques employed in this paper could also be applied to obtain Stein 
operators for independent products of a number of standard distributions. 


Table 1 Stein operators for product distributions. X ~ PB(ai, bi ..., am^bm), 
Y ~ PG(ri,..., r n , A) and Z ~ PN(JV, a 2 ) are mutually independent. 


Product P 

Stein operator Apf(x) 

Order 

X 

B ai ,...,a m f{x) - xB ai +b 1 ,...,a m + b m f(x) 

m 

Y 

Br 1 ,...,r n f[x) - X"xf(x) 

n 

Z 

a 2 A N f(x) - xf{x) 

N 

XY 

Bai ,... ,a m Br-y .. ,r n «/*(*^) ^ ^Ba^-\-b\ f 

m + n 

XZ 

& Bai ,... ,a m -An Bai ,... ,a m f (pB) 

+ • . ,a m + 6 m +&i — 1,... ,a m +b m — 1 f (*^) 

2m + N 

YZ 

a 2 Br 1 ,...,r n A N Br 1 ,...,r n f(x) - \ 2 "xf(x) 

2 n + N 

XYZ 

& Bai . .,a m B ri . ,r n An B r 1 f (^) 

A xB a ^- ? ... ,a m + 6 m -^ai +i>i — 1,... )Ct m + b m — 1 f (*e) 

2m + 2 n + N 


It can be seen that the product beta and product gamma Stein operators 
reduce to the classical beta and gamma Stein operators when m = 1 and 
n = 1, respectively, as was so in the normal case. In Section 2.2.2, we see 
that for certain parameter values the Stein operators for the products XZ 
and XYZ can be simplified to differential operators of lower order. We give 
a precise criteria under which this occurs. 

In Proposition 2.3, we also obtain a operator for the generalised gamma 
distribution which leads to the following PGG(?r,.. ., r n , A, q) Stein opera¬ 
tor: 

AIpgg f(x) = B ri ^ r „f(x) - ( q\ q ) n x q f{x ). (1.10) 

Taking q = 1 in (1.10) yields the product gamma Stein operator Ay f(x). 
Taking rq = • • • = rjy = 1, A = (y/2 cr) _1 and q = 2 in (1.10) gives the fol¬ 
lowing Stein operator for the product of N independent half-normal random 
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variables (\Z\ where Z ~ PN(N, cr 2 )): 

•ApHn/ (x) = (7 2 Ti f(x) - x 2 f(x), 

where x takes values in the interval [0,oo). By allowing x to takes values in 
M, we obtain the following PN(IV, a 2 ) Stein operator 

A z f{x) = u-'7' A f(x) - x 2 f(x), 

which differs from the PN(iV, cr 2 ) operator (1.9). Although, making the 
changes of variables g(x) = xf(x) we have that g'(x) = xf'(x) + f{x), 
and so 

A N g(x) = x~ 1 Tq g{x) = T^f{x), 

from which we recover the Stein operator (1.9). 

The product distribution Stein operators that are obtained in this pa¬ 
per have a number of interesting properties which are discussed in Remark 
2.6. However, despite their elegance, it is in general difficult to solve the 
corresponding Stein equation and bound the appropriate derivatives of the 
solution; further discussion is given in Remark 2.10. 

The classical normal, beta and gamma Stein equations are first order 
linear differential equations, and one can obtain uniform bounds for their 
solutions via elementary calculations. Uniform bounds are available for the 
first four derivatives of the solution of the PN(2, cr 2 ) Stein equation (Gaunt 
[13]), and in Proposition 2.8 we show that the k -th derivative of the solution 
of the PG(ri, r 2 , A) Stein equation is uniformly bounded if the first k deriva¬ 
tives of the test function h are bounded. Although, for all other cases of 
product distribution Stein equations we do not have bounds for derivatives 
of the solution. 

However, in Section 3, we consider a novel application of the product beta- 
gamma-normal Stein operator. In Section 3.2, we use the operator to obtain 
a differential equation that the product beta-gamma-normal p.d.f. must sat¬ 
isfy. This allows us to ‘guess’ a formula for the density function, which is 
then easily verified to be the correct formula via Mellin transforms. This 
formula is new, and obtaining it directly using the inverse Mellin transform 
would have required some quite involved calculations. From our formula we 
are able to obtain an expression for the characteristic function of the product 
normal-beta-gamma distribution, as well as estimates for the tail behaviour 
of the distribution. 

1.3 Outline of the paper 

We begin Section 2 by establishing some properties of the operators A/v 
and We then obtain characterising equations for mixed products of 
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beta, gamma and central normal random variables (Propositions 2.3, 2.4 and 
2.5), which lead to the operators of Table 1. In Section 2.2.2, we see that for 
certain parameter values simpler operators can be obtained. In Section 2.3, 
we consider a Stein equation for the product of two independent gammas. 
We solve the equation and show that the k- th derivative of the solution 
is uniformly bounded if the first k derivatives of the test function h are 
bounded. 

In Section 3, we obtain formulas for the p.d.f. and characteristic function 
of the product beta-gamma-normal distribution, as well as an asymptotic 
formula for the tail behaviour of the distribution. We use the product beta- 
gamma-normal Stein operator to propose a candidate formula for the p.d.f. 
and then verify it using Mellin transforms. 

In Appendix A, we prove some results that were stated in the main text 
without proof. Finally, Appendix B lists some basic properties of the Meijer 
G-function and modified Bessel functions that are used in this paper. 

Notation. Throughout this paper, we shall let T denote the operator 
Tf(x) = xf'(x) and An will denote the operator A^f(x) = x~ l T N f{x) = 
^(T^- 1 f{x))- We also let T r denote the operator T r f(x ) = xf'(x ) + rf(x ) 
and let S ri) ... jT . n denote the operator B rij ... jrn f(x) = T rn ■ ■ ■ T ri f(x). We shall 
let C n (I ) be the space of functions on the interval I with n continuous 
derivatives, and Cjf(I) will denote the space of bounded functions on I with 
n continuous derivatives that are all bounded. 

2 Stein operators for products of normal, beta and gamma 
random variables 

2.1 Preliminary results 

We begin by presenting some useful properties of the operators An and 

B ri ,...,r n ■ 

Lemma 2.1. The operators An and B r 1 r „ have the following properties. 

(i) The operators T r andT s are commutative, that is, T r T s f(x ) = T s T r f(x ) 
for all f G G 2 ( R). Thus, for all f G G n (M), B ri ^ Tn f{x) = £ Ml) ,...,rv (n) /(a0; 
where a is a permutation of the set {1,2,... ,n}. 

(ii) For all f € C n+N { M), the operators An and B rlt ... t r n satisfy 

ANB ri ^_^ Tn f (x) = 5 ri _|_i i ... irn _|_iAjv/(^)- (2-1) 

Proof, (i) The first assertion follows since T r T s f(x ) = x 2 f"(x) + (1 + r + 
s)xf'(x) + rsf(x ) = T s T r f(x), and the second assertion now follows imme¬ 
diately. 
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(ii) As Ai = we have A 1 T r f(x) = xf"(x) + (r + l)f'(x) = T r+1 A 1 f(x). 
Thus, on recalling that A^f{x) = f(x)) and using the fact that the 

operators T r and T s are commutative, we have 

A N B ri _rJ(x ) = A 1 T 0 N ~ 1 T ri • • • T r J(x) 

= A 1 T ri ---T rn T 0 N ~ 1 f(x) 

= T ri+1 A 1 T r2 ---T rn T 0 N ~ 1 f(x) 

= T ri+1 -..T rn+1 A 1 T 0 N ~ 1 f(x) 

= B ri+ i ) ^ t r n+ iA]\rf(x), 


where an iteration was applied to obtain the penultimate equality. □ 


The following fundamental formulas (Luke [26], pp. 24-26) disentangle 
the iterated operators An and B r i,For / € (7 n (R), 


N , . 

■4jv/W = E {"}x‘- 1 /“ l (i), 

k =: 1 

(2.2) 

TL 

(®) = ^2c k)n X k f {k \x), 

k =0 

(2.3) 


where {)!} = ^ X^=o( —J are Stirling numbers of the second kind 
(Olver et al. [30], Chapter 26) and 

c k,m = 7, ^2 — ^-{j + qa) + q r i), (2.4) 

j=0 J ■ i =1 

for (a)j = a(a + 1) • • • (a + j - 1), (a) 0 = 1. 

Applying (2.1) and (2.3) gives that, for / G C m+n+N (M.), 

^ai,...,a m ^4AfS?ii,...,6n/( a ') = B ai —l t ... t a m — lBb 1 ,...,b n f ( x ) 

= x- 1 T^ v S 01 _i ) ... ) 0 m _iS 6 li ... i 6 n /(x) 

m+n+TV 

= ^ 4,m+n+iVZ fe ~ 1 / (A:) ( a; ), (2-5) 

fc=l 

where the c k) m+n+N can be computed using (2.4). 
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2.2 Stein operators 

With the preliminary results stated, we are now in a position to obtain 
Stein operators for mixed products of beta, gamma and central normal 
random variables, which give rise to the product distribution Stein oper¬ 
ators of Table 1. From here on we shall suppose that the random variables 
X ~ PB(ai, b\ ... , a mi b m ), Y ~ PG(ri,..., r n , A) and Z ~ PN(iV, a 2 ) are 
mutually independent. We shall also let Apf(x) be the operator for the 
product distribution P, as given in Table 1. 

2.2.1 General parameters We firstly consider the case of mixed products of 
beta, gamma and central normal random variables with general parameter 
values. In Section 2.2.2, we look at particular parameter values under which 
we can obtain some slightly simpler formulas for product distribution Stein 
operators. We begin by recalling the product normal Stein operator that 
was obtained by Gaunt [13]. 

Proposition 2.2. Suppose Z ~ PN(JV, a 2 ). Let f € C N ( R) be such that 
E\Zf(Z)\ <oo and E|Z fc-1 / (fc) (Z)| < oo, k = 1,..., N. Then 

E[Azf(Z)} = 0. (2.6) 

We now present Stein operators for the product beta and product gener¬ 
alised gamma distributions; taking q = 1 gives a product gamma distribution 
Stein operator. 

Proposition 2.3. Suppose G ~ PGG(ri,..., r n , A, q). Let f € C n (R) be 
such that E\G q f(G)\ < oo and K\G k f^ k \G)\ < oo, k = 0 where 
/(°) = /. Then 


E [B ri _ r J(G) - (gATGV(G)] = (2.7) 

Proof. We proceed by induction on n and begin by proving the base case 
n = l. The well-known characterisation of the gamma distribution, given in 
Luk [25], states that if U ~ T(r/q, A), then 

E[Uf , (U) + (r/q-XU)f(U)]=0 (2.8) 

for all differentiable functions / such that the expectation exists. Now, if 

V ~ GG(r, X,q), then V = (A 1 ~ q U) 1 / q . Making the change of variables 

V = (A 1 -^) 1 ^ in (2.8) leads to the following characterising equation for 
the GG(r, X,q) distribution: 

Wf'iV) + (r - q\ q V q )f(V)] = 0 
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for all differentiable functions / such that the expectation exists. This can 
be written as E[T r f (V)—qX g V q f (V)\ = 0, and so the result is true for n = 1. 

Let us now prove the inductive step. We begin by defining W n = rr=i ^ 
where V) ~ GG(ri,X,q) and the Vi are mutually independent. We observe 
that ( T p f)(ax ) = T p f a (x) where f a (x) = f(ax), and so (B pu ... p J)(ax) = 
B pu ... Pl f a (x). By induction assume that ( qX q ) n EW%g(W n ) = EB ri} __^ rn g(W n ) 
for all g € C n (M) for some n > 1. Then 

(qX q ) n+1 EW q n+1 f(W n+1 ) 

= (qX q ) n+1 E[V q +1 E[W q fv n+1 (W n ) | V n+1 ]] 

= qX q E[V q +1 E[B ri _ r J Vn+1 (W n ) | V n+ x]\ 

= qX q E[V q +1 (B ri _ r J)(W n V n+1 )} 

= qX q E[E[V q +l {B ri _ r J Wn ){V n+l ) | W n ]] 

= E[E[W n V n+l (B ri ^ r jy(W n V n+1 ) +r n+1 f(W n V n+l ) \ W n ]] 

= EB ri ^^ rn+1 f(W n +i). 

Thus, the result has been proved by induction on n. □ 

Proposition 2.4. Suppose X ~ PB(ai, b \,..., a m , b m ). Let f € C' m ((0,1)) 
be such that R\X k fW(X)\ < oo and E\X k+1 f( k \X)\ < oo, k = 0 ,..., m. 
Then 

E[Axf(X)] = 0. (2.9) 

Proof. The proof is similar to that of Proposition 2.3, and we proceed by 
induction on m. Let W m = n'=i Xi where X, ~ Beta(aj, 6j) and the X t are 
mutually independent. The base case of the induction m = 1 is the well- 
known characterisation (1.4) of the beta distribution. By induction assume 
that E W m B ai+bu _ jam+bm g(W m ) = EB ait ,„ <arn g{W m ) for all g € C m (R ) for 
some m > 1. Then 


EfPm+l-Ba 1 +6 1 ,...,a m+1 +& m+1 /(B / TO +l) 

= E[X TO+ iE[W m 5 ai+ b 1) ... i0rri+ 6 m T 0rrt+1+ b m+1 /x m+1 (W m ) | X m+1 ]] 

= E[X m _|_iE[B ai) ,,, )am T arri+1 _(_6 m+1 /x m+1 (Wm) I X m+1 }] 

= E[X m+ 1 (T am+ 1 +W iB a i,..., aro /)(W m X m+1 )] 

= E[E[X m+ i(T arn+1+ b rn+1 B ai: _^ arn fw rn )(X m+ i) | W m ]] 

= E[E[X m _|_lW m (f? al) ... )am /w m ) , (X m+ i) + a m +i/(W m X m _|_i) | W m ]] 

= ES 0l) ... j0m+1 /(W m +i), 

and so necessity has been proved by induction on m. □ 


imsart-bjps ver. 2011/11/15 file: ngb_bjps_2016.tex date: December 14, 2016 


Products of normal, beta and gamma random variables 


11 


We now use the above product beta, gamma and normal Stein operators 
to obtain Stein operators for mixed products of such random variables. 

Proposition 2.5. Let X ~ PB(ai, b \..., a m , b m ), Y ~ PG(ri,..., r n , A) 
and Z ~ PN(1V, a 2 ) be mutually independent. 

(i) Let f G C m+n (U+) be such thatE\(XY)i f^(XY)\ < oo, j = 0,..., m+ 
n, and E\(XY) k+1 fW(XY)\ < oo, k = 0,..., m. Then 

E [Axyf(XY)\ = 0. 

(%%) Let f G C 2m+N (R) be such that E \(XZy- 1 (XZ)\ < oo , j = 

1,_, 2m + N, and E|(X Z) k+1 f( k \XZ)\ < oo, k = 0,..., 2m. Then 

E[Axzf(XZ)] = 0. (2.10) 

(in) Let f G C 2n+N (M.) be such that K\YZf(YZ)\ < oo and additionally 
E|(y Zf- 1 / (fc) (y z)\ <oo, k = l,...,2n + N. Then 

E[A Y zf(YZ)] = 0. (2.11) 

(iv) Let f G c 2m+2n+7V (R) be such that E {(XYZy^f^iXYZ)] < oo, 
j = l,...,2m + 2n + N, andE\{XY Z) k+l f^ k \XY Z)\ < oo, k = 0,... , 2m. 
Then 

E [A X Yzf(XYZ)\ = 0. (2.12) 

Proof. In our proof, we use the Stein operators of the product normal, prod¬ 
uct gamma and product beta distributions that were given in Propositions 
2.2, 2.3 and 2.4, respectively. We consider the four assertions separately. 

(i) Recall that ( T p f)(ax ) = T p f a (x) where f a (x) = f{ax), and therefore 
(Bpi,...pif)(ax) = Bp lt '" Pl f a {x). From (2.9) and (2.7) (with q = l),we now 
have 

A"E [XYB ai+bl _ am+bm f(XY)\ = \ n E[YE[XB ai+bl _ am+bm f Y {X) \ Y]] 

= A n E [YE[B ai ,...,a m f Y (X) | Y]] 

= A"E [YB ai _ a J{XY)\ 

= A"E[E [YB ai _ am f x (Y) | X]] 

= E[E[B ri _ rn B ai _ am f x (Y)\X}} 

= E [B ru ..., rn B au ... >am f(XY)], 


as required. 
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(ii) We begin by noting that, since Ajyf(x) = (Tq X f(x)), we have 
(An) f (ax) = aANf a (x). So from, (2.9) and (2.6), 

K[X ZB ai ^b 1 ,..., a m+ b mB ai +b 1 -l,...,a rn +b rn -lf (XZ)\ 

= E[ZE[XB ai+bu _ >am+bm B ai+bl ^ lt _ )am+bm _ 1 fz(X) | X]] 

= ^[ZE[B ai! '.. :am B ai+bl _ li ... tarn+brn _ l fz(X) | X]] 

= E[E[Z J B ai ,... !am S ai+fel _ 1) ... !am+6m _ 1 /. Y (Z) | X]] 

= a 2 E[E[XA N B air „ }arn B ai+bl _ 1} ___ tarn+brn _ 1 fx(Z) | X]] 

= (7 2 E[Xy4ivS ai| ... >0m S ai+&1 _i t ... )am+&m _i/(XZ)]. 

From Lemma 2.1 we can obtain that A N B ai ,... t am B a 1 +b 1 -i,...,am+b m -i = 
B ai +bi,...,a m +bm A NB au ... tam . Applying this formula and (2.9) yields 

E[XZB ai+bl ^__ jamYbrri B ai+bl _i ! ..., am +6 m _i/ (XZ)\ 

= cr 2 E[XB ai+bu ^^ arn+b7n A N B au ... yarn f(XZ)} 

= a^E[E[XB ai+bl ^_^ arn+brri ANB ait ,,, }arn fz(X) | Z]] 

= a 2 E[E[B ai _ am A N B ai _ arn f z (X) | X]] 

= a 2 E[B ai _ am A N B ai _ am f(XZ)\, 

as required. 

(iii) By a similar argument, 

A 2n E[YX/(XX)] = \ 2n E[ZE[Yfz(Y) \ X]] 

= X n E[ZE[B ri _ r J z (Y) | X]] 

= A"E[E [ZB ri _ r J Y (Z) | Y]] 

= a 2 X n E[E[YA N B ri _ r J Y (Z) \ Y]] 

= a 2 X n E[E[YA N B ri _ r J z (Y) | X]] 

= a 2 E[E(B ri _ rn A N B ri _ r J z (Y) | X]] 

= a 2 E[B ri ^ rn A N B ri ^ rn f(YZ)\. 

(iv) Applying (2.9) and (2.12) gives 

X 2n E[XYZB ai+bl ,...,a rn +bru B a 1 +b 1 -l,...,a rn +bm-lf(XYZ)] 

= A 2n E [XXE [XB ai+blt ___ tClim+bm B ai +bl -i, ...,a m +b m -ifyz(X) | XX]] 

= A 2 -E[YXE[ J B air .., aro S ai+6l _ 1 ,... )am+6m _ 1 /^(X) | XX]] 

= A 2n E[E[XX5 01i ... iara B 0l+fcl _ li ... >am+ftm _ 1 /x(XX) | X]] 

= a E[E[Xi? 7 . li ... irn Ajv-Bri,...,r„^ai,...,a m -Bai+6i-i,...,a m +6 m -i/x(F X) | X]] 

= a E[X5 ri| ... :rn Ajv^r 1) ..., rrl ^ ai ,...,a m -Ba 1 +6 1 —l,...,o m +6 m —l/(XXX)]. 
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We now interchange the order of the operators using part (ii) of Lemma 2.1 
and then use (2.9) to obtain 

A 2 "E [lFZB ai+6l ,..., am+6m B 0l+6l _ lr .., am+6ra _ 1 /(JF2)] 

= <J 2 K[XB ai+ i }1 ^^ am+ i )m B rit ." trn ANB ri ,..' trn B aij ... jam f(XYZ)\ 

= a E[E[X5 ai +b 1) _ iam +b m .B ri) ... irn ^7V-Bri,...,r n -Ba 1 ,...,a m ,/yz(^0 | YZJ\ 

= & M\K[B aij ' j(lm B rij '.. trn A nB ri ,... jrn B aij _ arn fYz{X) \ YZ]\ 

= U“E[S ai ,...,a m ^r 1 ,...,r n ^4Ar^r 1 ,...,r rl ^a 1 ,...,a m /(^h Z)\. 

This completes the proof. □ 

Remark 2.6. We could have obtained first order Stein operators for the 
product normal, beta and gamma distributions using the density approach 
of Stein et al. [41] (see also Ley et al. [24] for an extension of the scope 
of the density method). However, this approach would lead to complicated 
operators involving Meijer G-functions. We would expect that this would 
lead to various problems. Firstly, bounding the derivatives of the solution 
could still be a challenging problem, and these derivatives might not even 
be bounded. Moreover, the coupling techniques in the existing Stein’s method 
literature seem to be most effective when the coefficients take a simple form. 
Our Stein equations, on the other hand, are amenable to the use of couplings. 
Indeed, Gaunt [13] used a generalised zero bias coupling in conjugation with 
the product normal Stein equation to prove product normal approximation 
theorems. 

From the formulas (2.2) and (2.3) for the operators Ajy and B r it 
follows that the product Stein operators of Table 1 are linear ordinary differ¬ 
ential operators with simple coefficients. As an example, the Stein operator 
for the product XY Z can be written as 

2m+2n-\-N 2m 

A X Yzf(x) = a 2 ^2 a k,2 m +2n+NX k ~ 1 f (k \x)-\ 2n ^2fi k) 2mX k+1 f (k) {x), 

k =1 k =0 

where the a k ,2m+2n+N and fi k ,2m can be computed using (2-4). 

As discussed in the Introduction, Stein operators of order greater than 
two are not common in the literature; however, our higher order product 
Stein operators seem to be natural generalisations of the classical normal, 
beta and gamma Stein operators to products. It is interesting to note that 
whilst the product beta, gamma and normal Stein operators are order m, n 
and N, respectively, the operator for their product is order 2m + 2n + N, 
whilst one might intuitively expect the order to be m + n + N. The formula 
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(3.1) of Theorem 3.1 below for the p.d.f. for the product XYZ sheds light on 
this, and is discussed further in Remark 3.2. In Section 2.2.2, we shall see 
that for certain parameter values one can obtain lower order Stein operators 
for the product XYZ. For example, the operator decreases by m when b\ = 
■■■ = b m = 1, and this can also be understood from (3.1) and properties 
of the Meijer G-function; this is also discussed in Remark 3.2. However, 
for general parameter values, we expect that a Stein operator for XYZ with 
polynomial coefficients will be of order 2m + 2n + N; again, see Remark 3.2. 

2.2.2 Reduced order Stein operators By Lemma 2.1, we can write the Stein 
operators for the products XZ and XYZ as 



Axzf{x) = cr 2 x 1 B i 


and 



Axrzf(x) = a 2 x 1 B i 


With this representation, we can write down a simple criterion under which 
we can obtain Stein operators for the products XZ and XYZ with orders 
less than 2m + N and 2m + 2n + N respectively. For simplicity, we only 
consider the case of the product XYZ\ we can treat the operator for product 
XZ similarly. 

Define sets R and S by 

R — T b \,..., a m T b m , q>i F b\ 1,..., a m T b m 1}, 

iS — , ■ ■ •, flm, !)•••) Qm 1) Di • • • > r n , i 'i 1,..., r n 1,0,..., 0}, 

where it is understood that there are N zeros in S. Then if \R n S\ = t, the 
Stein operator AxYzffx) can be reduced to one of order 2m + 2n + N — t. 
To illustrate this criterion, we consider some particular parameter values, 
(i) &i = ■ ■ ■ = b m = 1: X is product of m independent U( 0,1) random 
variables when also a\ = ■ ■ ■ = a m = l. Here the Stein operator is 



fix) 

(2.13) 


where we used the fact that the operators T r and T s are commutative. Taking 
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g(x) = -B<n f (x) then gives the (m + 2n + lV)-th order Stein operator 



Ag{x) = (j 2 x 1 B, 


(2.14) 


It should be noted that the Stein operator (2.14) acts on a different class of 
functions to (2.13). To stress this point, we recall from part (iv) of Proposi¬ 
tion 2.5 that the operator (2.13) acts on all functions / € C 2rn+ 2n+Ar (M) 
such that 1 &\{XYZy~ 1 f^\XYZ)\ < oo, j = 1, ■ ■ ■, 2m + 2 n + N, and 
E\(XYZ) k+1 f<< k \XYZ)\ < oo, k = 0, ...,2m. Whereas, (2.14) acts on 
all functions g € C m+2n+iV (M) such that K\(XYZy~ 1 g^\XYZ)\ < oo, 
j = 1,..., m + 2n + N, and E| (XYZ) k+1 g^ (XYZ)\ < oo, k = 0,..., m. 

In the subsequent examples, we shall not write down the resulting lower 
order Stein operators, although they can be obtained easily by similar cal¬ 
culations. 

(ii) a\ + b\ = ■ ■ ■ = a m + b m = 1: X is a product of m independent arcsine 
random variables when also a\ = • • • = a m = 1/2. A Stein operator of order 
m + 2n + N can again be obtained. 

(iii) m = n = N, a\ + b\ = ■ ■ ■ = a m + b m = 1 and r\ = ■ ■ ■ = r n = 1, 
so that X and Y are products of m arcsine and Exponential 1) random 
variables respectively. A Stein operator of order 3m can again be obtained. 

(iv) m = n = N, a\ + b\ = ■ ■ ■ = a m + b m = 1 and n = ■ ■ ■ = r n = 2. A 
Stein operator of order 3m can be obtained. 

2.3 A Stein equation for the product of two gammas 

In general, for the product distribution Stein equations that are obtained in 
this paper, it is difficult to solve the equation and bound the appropriate 
derivatives of the solution. However, for the product normal Stein equation, 
Gaunt [13] obtained uniform bounds for the first four derivatives of the 
solution in the case N = 2. Here we show that, for the PG(ri,r 2 ,A) Stein 
equation, under certain conditions on the test function h, all derivatives of 
the solution are uniformly bounded. With a more detailed analysis than 
the one carried out in this paper we could obtain explicit constants; this 
is discussed in Remark 2.9 below. In Remark 2.10 below, we discuss the 
difficulties of obtaining such estimates for more general product distribution 
Stein equations. 

Taking q = 1 in the characterisation of the product generalised gamma 
distribution given in Proposition 2.3 leads to the following Stein equation 
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for the PG(ri,r 2 ,A) distribution: 

x 2 f"[x ) + (1 + n + r 2 )xf'(x) + (nr 2 - A 2 x)f(x) = h(x) - PG^ i r2 /i, (2.15) 

where PG^ T2 h denotes E h(Y), for Y PG(ri,r 2 ,A). The two functions 
a;-( r i +r ' 2 )/ 2 /P ri -r 2 (2Ay / x) and x~^ ri+r2 ^ 2 f ri _ r2 \(2Xy/x) (the modified Bessel 
functions I u (x) and K v (x) are dehned in Appendix B) form a fundamental 
system of solutions to the homogeneous equation (this can readily be seen 
from (B.12)). Therefore, we can use the method of variation of parameters 
(see Collins [5] for an account of the method) to solve (2.15). The resulting 
solution is given in the following lemma and its derivatives are bounded in 
the next proposition. The proofs are given in Appendix A. 

Lemma 2.7. Suppose h : M + —» M is bounded and let h(x) = h(x) — 
PG r 1 ,r 2 h- Then the unique bounded solution f : R + —>• M to the Stein equa¬ 
tion (2.15) is given by 

/M = - [ < (n+ra)/2 -T,_-.i(2^)fe(t) dt 

+ [t^+^-'K ri _ r2 (2\Smt)At (2.16) 

= - jf ‘ ( ' 1+r2)/2 - 1 ^,-. 2 | (2AVt)fc(0 dt 

- : ' f^ +rM2 -'K^ r ,(2XV-t)h(t)dt (2.17) 

Proposition 2.8. Suppose h G C^(M + ) and let f denote the solution (2.16). 
Then there exist non-negative constants Cq^, Cpfc,..., Ck,k such that 


<C 0)0 fh\\ and ||/ (fc) || <C 0 , k fh\\ + Y,C j , k \\h {j) l k>l. (2.18) 

j =i 

Remark 2.9. The solution f can be bounded by 

\f(x)\ < 2\\h\\ x{ri l 2)/2 J* t( ri + r2 )/ 2 - 1 1 K ri -r 2 (2Av / x)7|r 1 - r2 | (2A-\/t) 

- I| ri _ r2 |(2Av / x)Ar n -r 2 ( 2 A\/t)| dt, 
useful for ‘small’ x, and 

i/mi < 2 iaii ( 2AVt> dt 

+ 2 hii%lfT^/“ t(n+r2,/2_lA 'n-'-»( 2A Vt)dt, 
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useful for ‘large’ x. In the proof of Lemma 2.7, we use asymptotic formulas 
for modified Bessel functions to show that the above expressions involving 
modified Bessel functions are bounded for all x > 0. A more detailed analysis 
(see Gaunt [16] for an analysis that yields bounds for similar expressions 
involving integrals of modified Bessel functions) would allow one to obtain 
an explicit bound, uniform in x, for these quantities, which would yield an 
explicit value for the constant Co,o- By examining the proof of Proposition 
2.8, we would then be able to determine explicit values for all Cj } k by a 
straightforward induction. However, since we do not use the product gamma 
Stein equation to prove any approximation results in this paper, we omit this 
analysis. 

Remark 2.10. For the PN(2,cr 2 ) and PG(ri,r 2 ,A) Stein equations, one 
can obtain a fundamental system of solutions to the homogeneous equation 
in terms of modified Bessel functions. These functions are well-understood, 
meaning that the problem of bounding the derivatives of the solution is rea¬ 
sonably tractable. However, for product distribution Stein equations in gen¬ 
eral, it is more challenging to bound the derivatives, because the Stein equa¬ 
tion is of higher order and a fundamental system for the homogeneous equa¬ 
tion is given in terms of less well-understood Meijer G-functions (this can 
be seen from (B.8)), which do not in general reduce to simpler functions. See 
Gaunt [13], Section 2.3.2 for a detailed discussion of this problem for the 
product normal case. Obtaining bounds for other product distribution Stein 
equations is left as an interesting open problem, which if solved would mean 
that the Stein equations of this paper could be utilised to prove product, beta, 
gamma and normal approximation results. 


3 Distributional properties of products of beta, gamma and 
normal random variables 


3.1 Distributional theory 

Much of this section is devoted to proving Theorem 3.1 below which gives 
a formula for the p.d.f. of the product beta-gamma-normal distribution. 
Throughout this section we shall suppose that the random variables X ~ 
PB(ai, b \..., a m , b m ), Y ~ PG(ri ,... ,r n , A) and Z ~ PN(1V, a 2 ) are mutu¬ 
ally independent, and denote their product by W = XYZ. 
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Theorem 3.1. The p.d.f. ofW is given by 


p(x ) = KG 


2ra+2n+iV,0 

2m,2m-\-2n+N 


X 2n x 2 

2 2n+ N(T 2 


‘ ‘ ‘ n 

2 ) ’ 


ai+bi 
2 ' 
a l a m 

2 , . . . , 2 , 

ai+6i —1 

2 

In n -1 
•52’ 2 ’ ' 


fl|Tt +6 m 


Q-l 1 1 

2 J * * * ’ 2 ’ 

Q'm. J rbm. 1 
• ’ 2 

r '^" 1 0 0 ) ’ 
• • ? 2 ’ u ’ * * * ’ u / 


(3.1) 


where 

j- A" fV r(g, + 6 t ) -A- 2 r J 

V 2 2n+N/2 7T (n+N)/2 a Ll 2 b iT (o*) HrfoV 

i=l j = l 

We prove this theorem in Section 3.3 by verifying that the Mellin trans¬ 
form of the product XYZ is the same as the Mellin transform of the density 
(3.1). However, a constructive proof using the Mellin inversion formula would 
require more involved calculations. In Section 3.2, we use the product beta- 
gamma-normal characterisation (Proposition 2.5, part (iv)) to motivate the 
formula (3.1) as a candidate for the density of the product W. As far as 
the author is aware, this is the first time a Stein characterisation has been 
applied to arrive at a new formula for the p.d.f. of a distribution. 

Before proving Theorem 3.1, we note some simple consequences. The 
product normal p.d.f. (1.6) is an obvious special case of the master for¬ 
mula (3.1), and by using properties of the Meijer G-function one can also 
obtain the product beta-gamma density (1.7). 

Remark 3.2. Let us now recall the sets R and S of Section 2.2.2: 

R — T bi ,..., a m T b m , o>i + b\ 1,..., a m T b m 1}, 

S = {ai,... , flmj cli !)•••) 1, n, • • •, r n , n 1,... ,r n 1,0,..., 0}, 

where there are N zeros in set S. By property (B.l) of the Meijer G-function, 
it follows that the order of the G-function in the density (3.1) decreases by 
t if \R(1 S\ = t. This is precisely the same condition under which the order 
of the Stein operator Axyzf{x) decreases by t. The reason for this becomes 
apparent in Section 3.2 when we note that the density (3.1) satisfies the 
differential equation A* XY zPi x ) = 0> where A* XY z an a djoint operator of 
Axyz with the same order. Hence, the order of the Stein operator decreases 
precisely when the degree of the G-function in the density (3.1) decreases. 

As an example of this simplification, taking b\ = ■ ■ ■ = b m = 1 in (3.1) 
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and simplifying using (B.2), gives the following expression for the density: 


p{x) = KG 


ra+2n+iV,0 

m,m-\-2n-\-N 


X 2n x 2 

2 2n + N a 2 


ai+1 
2 ’ 

ffll — 1 Qm — 1 ri 

2 ) ■ ■ • J 2 ’ 2 ’ ' ' ' 


gro+1 
' ' ’ 2 

rn r 1-1 
2 > 2 ’ 


rVi—1 


, o,..., o ; 5 


where K is the normalizing constant. It is instructive to compare this with 
Example (i) of Section 2.2.2, in which a ( m+2n+N)-th order Stein operator 
was obtained for this distribution. 

The connection between the differential equation A* XY z'P( x ) = 0 an d the 
Stein operator AxYzf(x) also suggests that, for general parameter values, a 
Stein operator for XYZ with polynomial coefficients will be a ( 2m+2n+N )- 
th order differential operator. This is because the density (3.1) is a Meijer 
G-function that, for general parameter values, satisfies the (2m + 2n + N)-th 
order G-function differential equation (B.8). We expect this to be the case 
unless the sets R and S share at least one element. 

Finally, we record two simple corollaries of Theorem 3.1: a formula for 
the characteristic function of W and tail estimates for its density. We note 
that these formulas are new, but that there is quite an extensive literature 
on tail asymptotics for product distributions; see Hashorva and Pakes [21] 
and references therein; see also Pitman and Racz [35] for a recent neat 
derivation of the tail asymptotics for the density of the product of a beta 
random variable and an independent gamma random variable. 

Corollary 3.3. The characteristic function of W is given by 


</>(t) = MG 


2m-\-2n+N—l,l 
2m+l,2m+2n+AT—1 


x 2r 


2 2n +N-2 a 2 t 2 


n +1 

2 , • • 


1 2lL Am ±b m + 1 

2 ?*•• ? 2 ’ 

fli+1 dm + 1 ai dm 

2 ’ * * * J 2 ’2 , '** , 2’ 


Q1+&1 

2 ’ ' ' ' 

r n +1 n 
’ 2 » 2 


dm +b n 


, . . . , 2 ’ 2’’’'’ 


where 


M = 


1 

7r (n+N-l)/2 


n 


r(Oi + bi) -r-r 2 rj 1 

2 bi T(ai) 


Proof. Since the distribution of W is symmetric about the origin, it follows 
that the characteristic function (f{f) is given by 


= E[e itw ] 


E[cos(tW)] = 2 


f 

Jo 


cos (tx)p(x) dx. 
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Evaluating the integral using (B.5) gives 


— MC 2m+2n+A r ,l [ _”_ 

1\1 ^2m+2,2m+2n+N 1 2 2n + N ~ 2 a 2 t 2 


(. 


\ 2r ' 


1 ai +&i Q"m ±b 

m 

2 1 2 1 * • ' ' 2 5 
Qj Q-l 1 ftm 1 

2 ’ ' • ' ’ 2 ’ 2 ’•••’ 2 5 

ai+^l — 1 Qm+fom —1 fj 

2 ? ‘ * ‘ ’ 2 ’ u 


II 

2’***’ 2’ 2 ’•••’ 


where 


2iP0F _ 1 r(oi + bi) A 2 r i" 1 

|f| - nin+N -!)/2 2 6 *r(o i ) r(rj) ’ 2 n + JV/2 - 1 <r|t| ’ 


and simplifying the above expression using (B.2) and then (B.l) completes 
the proof. □ 


Corollary 3.4. The density (3.1) of the random variable W satisfies the 
asymptotic formula 


/ \ 2n x 2 \ l/(2n+iV) 'i 

p(x) ~ jV[x| a exp t - (2re + N) ( 22n+jVo . 2 j as |x| -> oo, 

where 

( 27r )(2n+JV-l)/2 / A 2n \ a/2 

^ “ (2n + AT) 1 / 2 ( V 2 2 «+ Ar a 2 J ^ 

with K defined as in Theorem 3.1, and 


a = 


2n + N 


1 — 3 n + N 


+E- 

l=i 


X>> 

l=i 


Proof. Apply the asymptotic formula (B.3) to the density (3.1). 


□ 


3.2 Discovery of Theorem 3.1 via the Stein characterisation 

Here we motivate the formula (3.1) for the density p of the product ran¬ 
dom variable W. We do so by using the product beta-gamma-normal Stein 
characterisation to find a differential equation satisfied by p. 

By part (iv) of Proposition 2.5 we have that 


IE[<7 B air ..,a rn Br 1 ,...,r n A]^B ri: ___ Xrl -B olj ... jam /(W) 

- X 2n WB ai+blt _ jarn+brn B ai+bl _ lj __ 4tarn+bm _ 1 f(W)] = 0 


(3.2) 
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for all / € C 2m+2n+iV (M) such that E|W fc_ 1 / (fc) (W)| < oo for 1 < k < 
2m+2n+N, and E\W k+1 f^\W)\ < oo for 0 < k < 2m. By using part (ii) of 
Lemma 2.1 and that A N f(x ) = and ^T k f{x) = T? +1 f(x), 


we can write 



B, 


T±,...,Vn j •••jQ'Ttt, 


On substituting into (3.2), we see that the density p(x) of W satisfies the 
equation 




CL1 ,...,CLm 


- A 2 n xB ai+bli ... iam+bm B ai+bl _i i ... jam+bm _i/(x)}p(x) dx = 0 (3.3) 

for all functions / in the class C p , which is defined by 

(i) / € C 2 m+ 2 n+iV (M); 

(ii) E|VL fc- 1 / (fc) (T / L)| < oo for 1 < k < 2m+2n+N and E\W k+1 f( k \\V)\ < 
oo for 0 < k < 2 m; 

(iii) x* +J+ 2 p(*)(x)/^^(x) —>■ 0 as x —> ±oo for all i,j such that 0 < i + j < 
2m — 1 ; 

(iv) x l+1 {x) f^Xx) —>• 0 as x —>• ±oo for all i,j such that 0 < i + j < 
2m + 2n + N — 1. 

It will later become apparent as to why it is helpful to have the additional 
conditions (iii) and (iv). Note that C p contains the set of all functions on R 
with compact support that are 2m + 2n + N times differentiable. 

We now note the following integration by parts formula. Let 7 € R and 
suppose that cj) and ^ are differentiable. Then 



— OO 




x 7+1 r cj)(x)-^-(x r il)(x)) dx 



(3.4) 


provided the integrals exist. 
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We now return to equation (3.3) and use (3.4) to obtain a differential 
equation that is satisfied by p(x). Using (3.4) we obtain 


f 

J —( 


0 

X p(x)-B ai _|_bi,...,a m _i+6 m _i-E>ai+&i-l,...,a m +6 m -l/ (x) 


L 

l 


xT 2 -a m -b m p(x)B ai+bl - lt _ tam+bm - 1 f(x)dx 
xT 2 - 0m - brn p(x)B ai+bl _ l ^ ani+bTn _ l f{x) dx, 


where we used condition (iii) to obtain the last equality. By a repeated 
application of integration by parts, using formula (3.4) and condition (iii), 
we arrive at 


xp(x)B a 1+bl ,... 

,am.+bm^ai+bi — l,...,a m +bm-l f(x) dx 
x f ( x )B3-a 1 -bi,...,3-a rn -b rn B2-ai-bi,...,2-a rn -b rn P( x ) dx. 

I 

By a similar argument, using (3.4) and condition (iv), we obtain 


J — C 


L 


JXJ 

P( x )B ai ,...,a rn Br 1 ,...,r n Br 1 +l,...,r n +lBai+l,...,a rn +lTi f (x) dx 
OO 

/ OO 

-OO 


X dx ^— a i’—>— a m^—ri,—,—r n Bi—r 1 ,...,l-r n Bi— ai: ___ i i- arn p(x)) dx. 

Putting this together we have that 

/ OO 

{(“ 1) a x T 0 i?- ai ,...,— am B- ri} ...,— r . n Si_ ri) ... ) i_ rn 5i_ air ,. i i_ am p(x) 

-OO 

- A 2n x5 3 _ ai _ 6li ... i3 _ am _ 6m S 2 _ ai _ felv .. i2 - am _ fem p(x)}/(x) dx = 0 


for all / € C p . Since the integral in the above display is equal to zero for 
all / € C p , it follows by a slight variation of the fundamental lemma of the 
calculus of variations (here we have restrictions on the growth of /(x) in the 
limits x —> ±oo) that p(x) satisfies the differential equation 

T 0 -B-ai,...,-a m -B-ri,...,-r n -Bl-ri,...,l-r Jl -E>l-ai,...,l-a m P(£) 

- (-l) Ar ^" 2 A 2n x 2 5 3 _ Q1 _ f , l! ... i3 _ am _ 6m S 2 _ ai _ fel) ... i2 _ am _ fem p(x) = 0. (3.5) 
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We now make a change of variables to transform this differential equation 
to a Meijer G-function differential equation (see (B.8)). To this end, let 
U = 2 ^+av • Then, x-^ = 2and p(y) satisfies the differential equation 


Tr\ B fo _ a ni B rjr^Bl—r-t l — r n Bl—a^ l — amPiy') 

2 5 ■ * * ? 2 2 2 —2—j***j—2— —2—)•••? v ' 

~ ( — 1) yB ‘3-a 1 -b 1 3-nm-tni B 2-a 1 -b 1 2 -am-bm Viv) = 0- (3-6) 

2 !•••> a 2 ’•••> 2 


From (B.8) it follows that a solution to (3.6) is 


Q.1+&1 Q-m +b m 

2 i • • • > 2 ’ 

QT Q-m 

2 , • • • , 2 5 

ai+fcl —1 Hm+tm-l 

2 ’ " ‘ ’ 2 

am —1 II III n —1 r n -l r, 

• i 2 > 2 ’ ‘ ‘ ‘ ’ 2 ’ 2 ’ ■ ■ ■ > 2 



n( v \ — fbf~<2rn+2n+N,0 ( 

P{y) ~ ULr 2m,2m+2n+N\ V 


a\ — l 
2 ’ * 


where C is a constant. Therefore, on changing variables, a solution to (3.5) 
is given by 


r) ( T ) — (b(-i2m+2n+Nfl I - - " 

p\X) 1 '- r 2m,2m+2n+N 1 2n+N a 2 


( X 2n x 2 


a-i+bi Urn ±b m 

2 i ■ ■ ■ > 2 ’ 

0-1 Om ai 1 dm 1 

2 i • • • > 2 ’ 2 •>■■■’> 2 > 

ai+fcl-1 Qm+fem-l 


II 
2 ’ 


In f'i-1 
’ 2 ’ 2 ’ 


r 2 , 
' n L 


, 0,...,0 


where C is an arbitrary constant. We can use the integration formula (B.6) 
to determine a value of C such that / M p(x)dx = 1. With this choice of C , 
p{x) > 0 and so p is a density function. However, there are 2m + 2n + N 
linearly independent solutions to (3.5) and whilst our solution p is indeed a 
density function, a more detailed analysis would be required to rigorously 
prove that it is indeed the density function of the product beta-gamma- 
normal distribution. Since a simple proof that p is indeed the density func¬ 
tion is now available to us via Mellin transforms, we decide to omit such an 
analysis. 


3.3 Proof of Theorem 3.1 

Firstly, we define the Mellin transform and state some properties that will 
be useful to us. The Mellin transform of a non-negative random variable JJ 
with density p is given by 

POO 

Mu(s) = EU s ~ 1 = / x s ~ 1 p(x)dx, 

Jo 
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for all s such that the expectation exists. If the random variable U has 
density p that is symmetric about the origin then we can define the Mellin 
transform of U by 

roo 

Mjj(s) = 2 / x s ~ 1 p(x) dx. 

Jo 

The Mellin transform is useful in determining the distribution of products 
of independent random variables due to the property that if the random 
variables U and V are independent then 

M uv (s ) = Mu{s)M v {s). (3.7) 

Proof of Theorem 3.1. It was shown by Springer and Thompson [38] that 
the Mellin transforms of X, Y and Z are 


M x (s) 

My(s) 


n 

3 = 1 


T (CLj -)- 6j) T(cij — 1 T s) 
T(cij) r(cjj T bj — 1 T s) 


1 




n 

3 =1 


T(rj - 1 + s) 

r(b) : 


A «s) = ^I_ 2 w( '- 1 )/ 2 ( r' _ 1 [r(|)]' v . 

Then, as the random variables are independent, it follows from (3.7) that 


Mxyz(s) = J 


j =i 


7T 


T (dj + bj) T(cij — 1 + s) 1 

T(oj) T(a_j + bj — 1 + s) 

L_ 2 K(,-l)/2 (T ,-l [r( , )] JV 


n 

3 = 1 


T(rj - 1 + s) 


r (rj 


(3.8) 


Now, let IT be a random variable with density (3.1). Since the density of 
W is symmetric about the origin, we have 

/*oo 

Mw(s) = 2 / x s ~ 1 p(x)dx 

Jo 
\ n 


_TT r («j + bj) tt 2 T i 

2 ‘ 2 n+N/ 2 7r (n+N)/ 2 a 11 2 b ^T(aP r(r,) 

7=1 v J 7=1 J 


n 


r(Oi+£)r(a L _L+£) 

W 


nr 


Tj + S 


2 n+ N /2 a 

A" 


r j ~ 1 


X 


[r(D] 


N 


(3.9) 


ii j-j. 

where we used (B.6) to compute the integral. On applying the duplication 
formula T(f )T(| + 1) = 2 1 ~ x ^/ttT(x) to (3.9) we can deduce that the ex¬ 
pressions (3.8) and (3.9) are equal. Hence, the Mellin transforms of W and 
XYZ are equal and therefore W and XYZ are equal in distribution. □ 
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Appendix A: Further proofs 


Proof of Lemma 2. 7. We begin by proving that there is at most one bounded 
solution to the PG(ri,r 2 ,A) Stein equation (2.15). Suppose u and v are 
bounded solutions to (2.15). Define w = u — v. Then w is bounded and is a 
solution to the homogeneous equation 

x 2 w"(x ) + (1 + r\ + r 2 )xw'(x ) + (nr 2 — A 2 x)w(x) = 0. (A.l) 

We now obtain the general solution to (A.l). We begin by noting that the 
general solution to the homogeneous equation 

x 2 s"(x ) + xs'(x ) — ( x 2 + (ri — r 2 ) 2 )s(x) = 0 


is given by s(x) = CH ri _ r 2 (i) + DI\ ri _ r2 \(x) (see (B.12)). Here, we have 
used the fact that K_ v (x) = K u (x) for any vEM and all x > 0, which can 
be seen immediately from (B.9). Thus, K\ ri _ r2 \(x) = A ri _ r 2 (x). A simple 
change of variables now gives that t(x) = CK ri _ r2 (2\y/x) + DI\ ri _ r2 \(2\y/x) 
is the general solution to 


x 2 t"(x) + xt'(x ) — (X 2 x + (n — r 2 ) 2 /4ft(x) = 0. (A.2) 


Substituting t(x) = x^ ri+r2 ^ 2 w(x) into (A. 2 ) now shows that w satisfies the 
differential equation (A.l), and we have that the general solution to (A.l) 
is given by 

w(x) = Aw\(x) + Bw 2 {x), 


where 


w\(x) = x ( r ’ 1 + r 2 )/ 2 ^ ri _ r , 2 ( 2\y/x) and W 2 {x) = x ^' 1+r2 ^ 2 I\ ri _ r2 \(2\y/x). 

From the asymptotic formulas for modified Bessel functions (B.10) and 
(B.ll), it follows that in order to have a bounded solution we must take 
A = B = 0, and thus w = 0 and so there is at most one bounded solution 
to (2.15). 

Since (2.15) is an inhomogeneous linear ordinary differential equation, we 
can use the method of variation of parameters (see Collins [5] for an account 
of the method) to write down the general solution of (2.15): 


. . f x W 2 (t)h(t) . . f x wi(t)h(t) . /an 

f(x) = -»,(x) f dt + w 2 (x) f dt, (A.3) 


where a and b are arbitrary constants and W(t) = W(wi,W 2 ) = wiw' 2 —W 2 W , 1 
is the Wronskian. From the formula W(K u (x), I u (x)) = x~ 1 (Olver et al. 
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[30], formula 10.28.2) and a simple computation we have W(wi(x),W 2 (x)) = 
\x~ 1 ~ ri ~ r " 2 . Substituting the relevant quantities into (A.3) and taking a = 
6 = 0 yields the solution (2.16). That the solutions (2.16) and (2.17) are equal 
follows because t^ ri ~ r2 ^ 2 ~ 1 AT n _ r2 (2A\/t) is proportional to the PG(ri, 7"2, X) 
density function. 

Finally, we show that the solution (2.16) is bounded if h is bounded. If 
r i 7^ f2, then it follows from the asymptotic formulas for modified Bessel 
functions (see Appendix B.2.3) that the solution is bounded (here we check 
that the solution is bounded as x | 0 using (2.16), and to verify that it is 
bounded as x —>• oo we use (2.17)). If ?q = r 2 , the same argument confirms 
that the solution is bounded as x —>• oo. To deal with the limit x | 0, we 
use the asymptotic formulas Iq(x) ~ 1 and Kq(x) ~ — log(x), as x 0, to 
obtain 




- I 0 (2X^)K 0 (2\Vt)} h(t) d t 


= a (ri+r a )/2 l t(n+r2)/2 “ 1 [ lQ g( X ) - 1 o g(*)] W) dt 

- Il ^ 11 iTo x(nlr 2 )/2 J Q t )/2 “ 1 [ lo g(x) - log(t)] d t 

= \\h\\ lim -- 1 = 4|1 ^ 11 

xj.0 ((n + r 2 )/2) 2 in + r 2 ) 2 


Therefore the solution is bounded when h is bounded. This completes the 
proof. □ 

Proof of Proposition 2.8. In this proof, we make use of an iterative approach 
that first appeared Dobler [8], and was then developed further in Dobler 
et al. [9]. Denote the Stein operator for the PG(ri,r 2 ,A) distribution by 
A ri ^ 2 ,\fi x )i so that the PG(r*i,r 2 ,A) Stein equation is given by 


A ri , r2 ,xfix) = hix). 


Now, from the Stein equation (2.15) and a straightforward induction on k, 
we have that 


x 2 f( k + 2 )( x ) _|_ ( ri + r 2 + 2k + 1 )xf^ k+ 1 \x) 

+ ((n + fc)(r 2 + k) — X 2 x)f( k \x) = h^ k \x) + kX 2 f^ k ~ 1 \x), 
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which can be written as 

Ar 1+ k,r 2 +k,\f ik \x) = h^ k \x) + k\ 2 f {k ~ l \x). 

Now, by Lemma 2.7, there exists a constant C ri)T2t \ such that 


— Cri,r2,A 


We also note that the test function h'(x)+\ 2 f(x) has mean zero with respect 
to the random variable Y ~ PG(ri + l, T 2 + 1, A), since by the product gamma 
characterisation of Proposition 2.3, 

E[h'(Y) + A 2 f(Y)} = E[A ri+k , r2+kjX f’(Y)\ = 0. 


With these facts we therefore have that 

ll/'ll < C n+hr 2 +ltX \\h'{x) + A 2 /(x)|| < C ri+hr 2 + 1 , x {\\h'\\ + A 2 ||/||) 

< Cri+l,r 2 +l,A(II^II + A^CV^^aIHI) • 

Repeating this procedure then yields the bound (2.18), as required. □ 

Appendix B: Properties of the Meijer G-function and 
modified Bessel functions 

Here we define the Meijer G-function and modified Bessel functions and 
state some of their properties that are relevant to this paper. For further 
properties of these functions see Luke [26] and Olver et al. [30]. 

B.l The Meijer G-function 

B.1.1 Definition The Meijer G-function is defined, for 2 € C \ {0}, by the 
contour integral: 

CL \,... , CLp \ 

bl, ■ ■ ■ ,b q ) 

j_ r c+io ° z _ s wjLi r ( g + b i) n?=i r (! - - g ) d 

2vr iJc-ioo Z Ilj=n +1 r (» + °j) rij= m +i r (! - -s) 

where c is a real constant defining a Bromwich path separating the poles of 
r(s + bj ) from those of T(1 — dj — s) and where we use the convention that 
the empty product is 1. 
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B.1.2 Basic properties The Meijer G-function is symmetric in the param¬ 
eters a±,..., a n ; a n+ i,..., a p ; b\, ..., 6 m ; and b m+ 1 ,... , b q . Thus, if one the 
a/s, j = n + 1,... ,p, is equal to one of the b^s, k = 1,..., m, the G-function 
reduces to one of lower order. For example, 


r»m,n 


Qj\ 5 • • • 5 &p— 1 ) &1 

b 1: ...,b q 


_ wm— l,n 

— ^p- 1 , 9-1 


The G-function satisfies the identity 


cs~im,n I 
* P,Q 


{, 

a i, • 

■ , a P ^ 

= G™’ n (z 

V 

bi,. 

■,bq) 

p,g y 


CLl ,... , dp— 1 

b 2 ,...,b q 


a\ + c,... ,a p + c 
bi + c,... ,b q + c 


m,p,q > 1 . 

(B.l) 

(B.2) 


B.l.3 Asymptotic expansion For x > 0, 

-x J exp ( — crx 1/,<T ), as x —>• oo, (B.3) 


G 9,0 

p,q 


X 


ai,...,a p \ _ ( 2 vr)( ,T 1)/2 ^ _ i/^ 


bi,---,b q 
where a = q — p and 


1/2 


6 > = - 
a 


a 


1 \ 1 — a 


v 

- V a r 

i=l i =1 


+ bi - 1 


B.1.4 Integration 


f 


e ux G 1 ax 


°1, ‘ ‘ ' ’ a P \ , 1 „. _ , —ls~<m,n+l t a 

p+1 ’’ 


UJ 


0 , oi,..., a p 
bi,--.,b q 


(B.4) 

For the conditions under which this formula holds see Luke [26], pp. 166-167. 

For a > 0, 7 > 0, aj < 1 for j = 1,..., n, and bj > — \ for j = 1,..., m, 
we have 


f 

■Jo 


cos('yx)G 7 T’ n ( ax 


_ fZ.. .-l/om.n+l /4a 

-V7T7 G p+ 2,q 


bi,...,b q 
2 , ® 1 , • • • 


dx 


& 1 , 


■ j 1 0 

! bq 


(B.5) 


The following formula follows from Luke [26], formula (1) of section 5.6 
and a change of variables: 


f 

■Jo 


x s ~ l G™: q n ( ax 2 


ai j•■ • j 


dx 


a 


bi,...,b q 

g n;7 =1 r(6, + f)n- =1 r(i-a,-f) 


2 nLm+i r(l - bj - f) IlLn+i r ( a l + I) 


(B. 6 ) 

(B.7) 


For the conditions under which this formula is valid see Luke, pp. 158-159. 
In particular, the formula is valid when n = 0, 1 < p+l < m < q and a > 0. 
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B.1.5 Differential equation The G-function f(z) = \ a b f"’ a b ') satis¬ 

fies the differential equation 

(-iy- m -n Z B 1 _ ai _ 1 _ ap f(z) - B_ bl __ bq f(z) = 0, (B.8) 

where B riy ^ rn f(z) = T rn ■ • • T ri f(z ) for T r f(z ) = zf'(z) + rf(z). 


B.2 Modified Bessel functions 

B.2.1 Definitions The modified Bessel function of the first kind of order 
n € M is defined, for all x € R, by 


= Y 


1 


k =0 


T(i/ + k + l)fc! 


x\ v + 2k 
2/ 


Ux) = r"G 


■ -^2,0/ x 


0,2 

a'4x) = 

B.2.3 Asymptotic expansions 
Iu(x) 


V V 

2’~2 y ’ 


x € 


1 


i/ ^ , 

2’“2>’ I>0 ' 


X | 0, 


Ky{x) 

Iy{x) 

Kjx) 


T(z2 T- 1) V 2 / 
r2i^i- 1 r(|^|)x~i^i, x|o,!//o, 

1 — log x, x | 0, i/ = 0, 

e* 

, X —>• oo, 


\/2vrx 
7r 


(B.9) 


The modified Bessel function of the second kind of order ^ £ R is defined, 
for x > 0, by 

/•OO 

K„(x)= Q— x c°sh(t) cosh(t't) dt. 

Jo 

B.2.2 Representation in terms of the Meijer G-function 

„2 


(B.10) 

(B.ll) 


B.2.4 Differential equation The modified Bessel differential equation is 

x 2 f"(x) + xf(x) - (x 2 + v 2 )f{x) = 0. (B.12) 

The general solution is /(x) = AI u (x) + BK v {x). 
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